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SUMMARY 


This  report  describes  the  basis  and  use  of 
the  Fortran  4  program  ADSYM  (Automobile  Defog/ 
defrost  SYstem  Model).  The  scope  and  limitations  of 
the  model  and  the  use  of  the  program  are  illustrated 
by  examples.  Possibilities  for  future  development 
are  briefly  discussed.  A  complete  listing  of  the  pro¬ 
gram  is  given  in  an  appendix. 
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ADSYM:  A  FORTRAN  MODEL  OF  AN  AUTOMOBILE 
DEFOG/DEFROST  SYSTEM 


INTRODUCTION 

An  adequate  field  of  vision  for  the  driver  is  evidently  a  basic  necessity  for 
the  safe  operation  of  an  automobile.  The  problem  of  clearing  and  maintaining  an 
adequate  field  of  vision  is  one  of  particular  importance  in  the  rather  severe  weather 
conditions  typical  of  a  Canadian  winter. 

A  series  of  tests  on  Canadian-built  cars  has  recently  been  conducted  for  the 
Ministry  of  Transport  by  the  Structures  and  Materials  Laboratory  of  the  National 
Aeronautical  Establishment  and  the  Low  Temperature  Laboratory  of  the  Division  of 
Mechanical  Engineering.  The  purpose  of  these  tests  was  to  provide  data  related  to 
CMVSS  103  on  windshield  defrosting  and'defogging  and  to  possible  future  requirements 
in  respect  of  rear  windows.  This  work,  and  some  related  earlier  experiments  are 
described  in  References  1  to  3. 

Cold  chamber  tests  of  this  type,  currently  required  by  S.A.E.  and  U.S. 
Federal  Government  Specifications  are,  however,  both  expensive  and  time  consuming. 
For  the  extensive  parametric  studies  involved  in  the  development  of  either  new  defrost 
and  defog  systems,  or  of  new  performance  standards,  faster  and  more  economical 
techniques  are  needed  to  reduce  the  time  required  in  the  cold  chamber.  The  value, 
in  this  context,  of  a  mathematical  model  of  an  automobile  defog/defrost  system  was 
recognized  by  Blatt  and  collaborators*4’ .  In  their  extensive  study  for  a  group  of  U.  S. 
Federal  Government  agencies,  they  developed  an  analytical  model  which  provides,  in 
principle,  a  basis  both  for  independent  parametric  studies  and  for  correlating  experi¬ 
mental  results  obtained  under  differing  conditions. 

The  practical  utility  of  the  original  model  was  rather  severely  limited  by  its 
analytical  form;  the  present  report  therefore  describes  a  version  of  the  model  pro¬ 
grammed  in  Fortran  and  known  as  ADSYM  (Automobile  Defog/defrost  SYstem  Model). 
While  the  programmed  version  adheres  closely  to  the  basic  equations  of  the  original 
model,  some  improvements  have  been  made  by  substituting  or  including  numerical 
methods  rather  than  analytical  ones. 

No  attempt  is  made  here  to  reproduce  the  extensive  discussion  and  analysis 
given  by  Blatt  et  al.  in  developing  the  model.  In  general,  the  basic  equations  are  merely 
restated  with  the  major  assumptions  relating  to  their  derivation  and  use.  Where 
ADSYM  differs  in  a  major  respect  from  the  original  model  however,  the  additional  or 
alternative  derivation  is  fully  stated. 


SYSTEM  MODEL 
General  Description 

The  complete  system,  comprising  the  environment,  the  vehicle  and  its  occu¬ 
pants,  is  represented  in  the  model  as  an  assembly  of  five  subsystems.  These  are: 

(i)  Vehicle  exterior  subsystem 
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^ii}  Engine  and  heater  subsystem 
(i;*'}  Vehicle  interior  subsystem 
(iv>  Interior  boundary  layer  subsystem 
(v)  Windshield/backlight  subsystem. 

Ambient  temperature  and  humidity  are  defined  by  the  vehicle  exterior  sub¬ 
system.  These  parameters  in  turn  define  the  occurrence  and  form  of  precipitation  and 
hence  the  heat  and  mass  transfer  mechanism  on  the  exterior  surface  of  the  windshield/ 
backlight. 


The  parameters  of  the  engine  and  heater  subsystem  define  the  transient 
temperature  behaviour  of  the  engine  coolant  and  hence  of  the  heater  exit  air. 

The  vehicle  interior  air  temperature  is  related  directly  to  the  heater  exit  air 
temperature  by  the  vehicle  interior  subsystem.  The  interior  moisture  concentration 
includes  the  effects  of  occupants'  breathing  and  vapour  generation,  when  appiopriate, 
from  wet  clothing. 

The  interior  boundary  layer  subsystem  may  be  specified  in  either  of  two 
forms.  The  first  represents  the  heat  and  mass  transfer  processes  which  occur  when 
a  hot  defroster  jet  mixes  with  the  cooler  and  wetter  interior  air.  The  second  form 
represents  free  convection  between  the  interior  air  and  the  surface  of  the  windshield/ 
backlight. 


A  windshicld/backlight  of  either  laminated  or  isotropic  construction  may  be 
specified  in  the  last  subsystem.  A  uniform  heating  rate  at  a  surface  or  interface  may 
also  be  specified,  to  represent  an  electrical  resistance  element.  The  subsystem  con¬ 
siders  the  heat  flow  through  the  windshield/backlight  under  the  time-dependent  boundary 
conditions  determined  by  the  first  four  subsystems.  The  surface  temperatures  at  a 
specified  location  are  determined  and,  with  the  surface  moisture  balances,  lead  to  a 
prediction  of  the  surface  states  at  that  location. 

Vehicle  Exterior  Subsystem 

The  vehicle  exterior  subsystem  defines  the  ambient  temperature  and  mois¬ 
ture  concentration  throughout  the  temperature  transient.  In  consequence,  it  also  de¬ 
fines  the  initial  surface  condition  and  heat  transfer  mechanism  for  the  exterior  of  the 
windshield/backlight.  The  table  below  shows  the  conditions  which  result,  in  the  present 
model,  from  specifying  the  various  possible  combinations  of  temperature  and  humidity. 


Ambient 

Temperature 

°F 

Relative 

Humidity 

% 

Initial  Surface 
Condition 

Heat  Transfer  Mechanism 

All 

<100 

Dry 

Forced  convection 

<25 

) 

.09  lb/ft2  of  ice 

Forced  convection 

25-32 

^100 

.  09  lb/ft2  of  ice 

Rainwater  film  cooling 

>32 

) 

wet 

Rainwater  film  cooling 
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For  icing  conditions,  the  surface  loading  assumed  is  that  specified  in 
SAE  J902a  for  windshield  defrosting  performance.  At  ambient  temperatures  less  than 
25°F,  the  effect  of  precipitation  on  the  heat  transfer  process  is  ignored.  A  precipita¬ 
tion  rate  of  one  inch  per  hour  has  been  assumed  in  rainfall. 

The  heat  transfer  coefficient  for  forced  convection  is  estimated  from  the 
expression  for  a  uniformly  heated  flat  plate  with  a  turbulent  boundary  layer: 


Nu  =  =  0.0296  Re8  Pr  33 

k 


0) 


where 


Nu 

is 

the 

h 

is 

the 

X* 

is 

the 

k 

is 

the 

Re 

is 

the 

Pr 

is 

the 

Utilizing  the  approximate  analogy  between  heat  and  mass  transfer  processes, 
the  corresponding  mass  transfer  coefficient  is  given  by: 


Sh  =  =  0.0296  Re8  Sc 33 

where  Sh  ‘.s  the  Sherwood  number 

g  is  the  local  mass  transfer  coefficient 
D  i*?  the  diffusion  coefficient  for  water  in  air 
P  is  the  density  of  air 
Sc  is  the  Schmidt  number  for  air 


(2) 


Combining  Equat’ons  1  and  2  the  relation 

_h_ 

CP 

is  obtained. 


g  = 


(10 


67 


(3) 


Considering  the  temperature  dependence  of  the  several  properties  of  air 
involved  in  Equations  1  and  3,  it  appears  that  the  heat  and  mass  transfer  coefficients 
may  be  assumed  independent  of  temperature  with  little  loss  of  accuracy  over  the  ambi¬ 
ent  range  of  interest.  Accordingly,  the  properties  of  dry  air  at  15°F,  laken  as 

k  =  0.135  x  10' '  Btu/ft  hr  °F 
v  =  0.135  x  10'3  ftVsec 
Pr  =  0.710 
cp  =  0.240  Btu/lb  °F 
D  =  0.234  x  10'3  ft2/sec 
have  been  used  for  computation 
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In  rain,  the  exterior  surface  of  the  windshield/backlight  is  assumed  to  be 
cooled  by  a  rainwal^i  film.  The  film  heat  transfer  coefficient  is: 

h  -  T  <4> 

where  k  is  the  thermal  conductivity  of  water  and  t  is  the  film  thickness. 

In  Appendix  A  a  simple  model  of  >.he  flow  over  the  windshiold/backlight  is 
developed  which  leads  to  a*'  estimate  of  the  film  thickness  t.  The  following  properties 
of  water  at  32°F  have  been  assumed  in  computation: 

k  =  0.345  Btu/ft  hr  **F 
"  =  0.1925  x  10'4  ft2/sec 
'  ^  62.57  lb/ft3 

The  v.  erature-dependence  of  the  properties  of  water  has  been  ignored  in  the  small 
temperature  range  of  practical  interest. 

The  details  of  the  heat  and  mass  transfer  processes  at  the  interface  between 
the  film  and  the  ambient  air  have  been  ignored.  The  effects  on  the  heat  transfer  coef¬ 
ficient  that  could  reasonably  be  represented  in  the  model  are  probably  quite  small; 
rainfall  is  the  dominant  mass  transfer  mechanism. 

Engine  and  Heater  Subsystem 

Consideration  of  the  heat  balance  for  the  engine  and  heater  subsystem  leads 
to  a  first  order  differencial  equation  for  the  coolant  temperature  Tc  at  time  t  with 
solution: 

(Tt-T„)  =  Tco  +  (Tt„-Tto)  (l-e‘,,T)  (5) 

where  Tco  =  (Tc-T0)!iQ 
Tc„  =  (Tc-T0)(_ 

Tc  =  ambient  temperature 
r  =  characteristic  time  of  subsystem 

Operation  of  the  thermostat  imposes  the  limit: 

T  <  T 

c  '  cmax 

where  TcriOJ(  is  the  coolant  temperature  at  which  the  thermostat  opens. 

The  temperature  of  the  defroster  jet  as  it  enters  the  car  is  given  by 

Tjet  ~T0  =  eH(Tc-T0)  (6) 


where  eH  is  the  heater  effectiveness,  assumed  constant. 


IThe  parameters  r,  Tc<>  and  eH  are  taken  to  be  independent  of  ambient  tem¬ 
perature.  They  are  determined  by  fitting  curves  of  the  general  form  of  Equat;on  5  to 
experimental  measurements  of  the  coolant  and  defroster  jet  temperatures.  The  values 
£  of  r,  Tt„  and  eH  are  dependent  on  other  parameters,  external  to  the  model,  such  as 

I  engine  idling  speed  and  heater  control  valve  position.  Consequently  they  are  defined 

or.ly  for  specified  sets  of  such  external  parameters. 

The  moisture  concentration  in  the  defroster  jet  is  assumed  to  be  the  same 
as  the  amo.ei.t  value,  since  no  moisture  is  added  to  the  air  passing  through  the  heater. 

Vehicle  Interior  Subsystem 


Consideration  of  the  interior  heat  balance  leads  to  an  approximate  expression 
for  the  interior  temperature  TCOf : 

TCQ,-T0  =  w*et(Tiwl-T0)  (7) 

where  w*  e,  is  a  dimensionless  constant  involving  the  heater  mass  flow,  leakage  flow 
into  the  vehicle  and  a  thermal  l^ss  parameter.  It  is  assumed  independent  of  ambient 
temperature  and  may  be  determined  for  a  particular  vehicle  from  measurements  of 
TCOI  and  T|et  during  a  temperature  transient. 


Consideration  of  the  moisture  balance  leads  to  an  approximate  expression 
for  the  interior  moisture  concentration  m,„,: 


N0  JO- 16  +  Gocc(mC0I  sot  -mD, 


mc0t-mo  = 


>] 


w  +  N  G 

jm  A  occ^occ 


(8) 


where 


is  the  ambient  moisture  concentration 
is  the  number  of  occupants 

0.16  is  the  vapour  generation  rate  in  lb/hr  due  to  one  person1  s  breathing 
is  the  transfer  coefficient  for  wet  clothing 

is  a  composite  mass  flow  determined  mainly  by  the  heater  mass  flow. 
Interior  Boundary  Layer  Subsystem 

The  tenmera.ure  and  moisture  concentration  profiles  in  a  defroster  jet  mix¬ 
ing  with  the  vehicle  interior  air  are  assumed  to  be  given  by  the  empirical  expression: 


mo 

is 

Ncee 

is 

0.16 

is 

Gocc 

is 

w.m 

is 

T*-Tc„,  mx-mCO( 


T,e,  ~T<or  m0-m« 


=  « 


\  1.6  )"2 

(1  +  0.15(x/z’|el)  | 


(9) 


is  the  maximum  temperature  in  the  jet  at  distance  x  from  the  exit  slot 
is  the  corresponding  maximum  vapour  mass  fraction 
is  the  decaj  function 

rz|e, 

ratio  of  exit  slot  length  to  width  of  windshield 
width  of  exit  slot  (normal  to  windshield). 
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The  local  heat  transfer  coefficient  h  is  obtained  from  the  related  empirical 
expression: 

hx*  65  33 

Nil  ~  =  c  Re  Pr  (101 

where  Nu  is  the  Nusselt  number 
x'  =  x  +  z'|t,  /0. 15 
k  is  the  thermal  conductivity  of  air 
c  =  0.405  z'|et325 

Re  is  the  Reynolds  number  based  on  x'  and  the  initial  jet  velocity  viel 
Pr  is  the  Prandtl  number  for  air. 

The  analogy  between  heat  and  mass  transfer  processes  again  leads  to 
Equation  3  for  the  mass  transfer  coefficient: 


where  in  this  case  h  is  given  by  Equation  10. 

Considering  the  temperature  dependence  of  the  properties  of  air  involved  in 
Equation  10,  it  appears  that  the  interior  heat  transfer  coefficient  is  more  nearly  inde¬ 
pendent  of  temperature  than  the  corresponding  exterior  coefficient  given  by  Equation  1. 
However,  while  the  ambient  temperature  remains  constant,  the  defroster  jet  temper¬ 
ature  varies  markedly  with  time  and  over  a  greater  range  than  that  of  interest  in 
specifying  ambient  conditions.  For  this  reason,  the  temperature  dependence  of  the 
conductivity  and  kinematic  viscosity  of  air  has  been  retained  in  this  part  of  the  model. 
Applying  Newton's  interpolation  formula  to  Eckert  and  Gross's  Table  A-4  of  Reference  5, 
the  following  expressions  are  obtained  for  the  temperature  range  -10  to  +170°F: 

v  =  0.1218  +  0.00052  (T+10)  -  0.0000004038  (T+10)  (T-80)  +  ... 

k  =  0.1287  +  0.000254  (T+10)  -  0.0o000006l72  (T+10)  (T-80)  +  ... 

where  T  is  the  air  temperature  in  degrees  Fahrenheit 
v  is  the  kinematic  viscosity  in  10'3  ft2/see  units 
k  is  the  conductivity  in  10 ' 1  Btu/ft  hr  °F  units. 

The  remaining  relevant  properties  are  assumed  to  be  constants  with  values 

at  80°F: 

cP  =  0.240  Btu/lb°F 
Pr  =  0.708 
Sc  =  0.577. 


Ir  the  absence  of  an  air  jet  blowing  over  the  interior  surface,  heat  and  mass 
transfer  between  the  surface  and  the  vehicle  interior  are  assumed  to  occur  by  free 
convection.  For  horizontal  or  vertical  plates  in  turbulent  free  convective  flow,  the 
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ASHRAE  Handbook  of  Fundamentals^’,  gives  the  relation: 

N«  =0.13  (Gr  Pr)1/3  (11) 

where  Nu  is  the  Nusselt  number 
Gr  is  the  Grashof  number 
Pr  is  the  Prandtl  number. 

At  temperatures  within  the  range  of  interest,  the  simplified  relations 

h  =  0.19  (To-T,)'*3  (vertical  plate) 

h  =  0.22  (Ta-Ts),'/3  (horizontal  plate) 

are  given  for  air.  For  most  angles  of  inclination  of  the  windshield/backlight  it  is 
assumed  to  be  sufficiently  accurate  to  take: 

h  =  0.205  (TtCf-T,)  33  (12) 

where  h  is  the  average  heat  transfer  coefficient 
Tcor  is  the  vehicle  interior  air  temperature 
T,  is  the  windshield/backlight  interior  surface  temperature. 

The  analogy  between  heat  and  mass  transfer  processes  allows  us  to  write: 

Sh  =  0.13  (Gr’  Sc) 1,3  (13) 

where  Sh  is  the  Sherwood  number 

Gr’  is  a  modified  Grashof  number  based  on  density  potential  rather  than 
temperature 

Sc  is  the  Schmidt  number 

Combining  Equations  11  and  13  a  relation  between  heat  and  mass  transfer  coefficients 
in  free  convection 


can  be  obtained  which  is  analogous  to  Equation  3  for  forced  convection.  Now 


Gr' 

Gr 


P  /P,-l 

cor  1 

0<Ttol-T,) 


(15) 


from  the  basic  definitions  of  the  two  Grashof  numbers'55,  where  the  numerator  is  the 
density  potential  between  the  surface  and  the  vehicle  interior  and  (i  is  the  volumetric 


, .  ^-,  ^ , : ,.-  ,f?^v:  jgSLiO^  i 
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thermal  expansion  coefficient  of  the  air.  For  small  mass  fractions  of  water  vapour, 
the  density  potential  may  be  expressed  in  terms  of  the  mass  fractions  as: 

P 

-  1  =  0.611  (mrmCQt)  +  0(m2)  (16) 

where  m,  is  the  saturation  mass  fraction  at  the  surface  and  mcor  is  the  mass  fraction 
of  water  vapour  in  the  vehicle  interior.  A  representative  value  of  /3  0.00195°F~'  is 

consistent  with  the  previous  development  of  Equation  12.  Taking  values 


cp  =  .240  Btu/lb°F 
Pr  =  .708 
Sc  =  .577 


and  substituting  also  from  Equations  15  and  16  into  14  we  obtain  finally: 

g  =  6.6  (m,-mcor) 33 


(17) 


Windshield/Backlight  Subsystem 


The  heat  balance  for  the  windshield/backlight,  neglecting  in- plane  conduction, 
is  expressed  by  the  one-dimensional  heat  conduction  equation  with  time-dependent 
boundary  conditions  at  the  interior  and  exterior  surfaces.  In  the  original  model,  this 
equation  was  solved  for  the  surface  temperatures  in  a  laminated  windshield,  using 
rather  crude  approximations  to  the  temperature  distribution  and  windshield  properties. 
In  the  present  model,  which  is  not  constrained  by  the  need  to  produce  analytical  rela¬ 
tions  for  the  variables  of  interest,  a  laminated  windshield/backlight  may  be  represented 
as  easily  as  one  which  is  isotropic.  Further,  a  uniform  interlaminar  heating  rate, 
representative  of  an  electrical  resistance  heating  element,  may  readily  be  specified. 

In  this  section  of  the  report,  the  basic  equations  are  stated  as  concisely  as  possible. 
In  Appendix  B,  the  equivalent  finite  difference  equations  are  given  and  the  method  of 
solution  is  described. 


The  equation  governing  heat  conduction  within  a  given  lamina  is: 


9T 

9t 


_k_  a^T 
cp  3z2 


(13) 


where 


T  is  the  local  temperature 
t  is  the  time  from  initial  thermal  equilibrium 
k  is  the  thermal  conductivity  of  the  lamina 
c  is  the  specific  heat  of  the  lamina 
p  is  the  density  cf  the  lamina 
z  is  the  normal-to-plane  co-ordinate 
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At  the  interface  between  adjacent  laminae,  continuity  of  heat  flow  requires  that: 


= 


where  Hfh  is  the  Interlaminar  heating  rate  and  the  subscripts  1  and  2  indicate  the  parti¬ 
cular  laminae  adjacent  to  the  interface.  Similarly,  at  the  two  surfaces,  the  boundary 
conditions  are  defined  by: 

8T,  .  \ 

k>  ft;  =  h-<T>-T-)  “  H>- 

\  (20) 
0Tn 

K„  9 —  =  he(Te-Tn)  +  H,e  / 

where  subscripts  1  and  n  indicate  respectively  the  innermost  and  outermost  laminae,  h, 
and  he  are  the  interior  and  exterior  heat  transfer  coefficients,  T,  and  Te  the  corres¬ 
ponding  temperatures  in  the  adjacent  air,  H, ,  and  H,  e  the  total  superficial  heating  rates 
from  all  other  sources. 

In  general: 

Hs  =  Hrh  +  Hmf  -  Hmelt  (21) 


Hrh  is  the  superficial  heating  rate  due  to  an  electrical  resistance  element,  a  known 
constant.  Hmt  is  the  superficial  heating  rate  due  to  mass  transfer  to  the  surface: 


Hmt.,  =  g,(m,-m,  Jot)  hv 

Hmt.e  —  6e(^e—^n,jot) 


(22) 


where  g|  and  g„  are  respectively  the  interior  and  exterior  mass  transfer  coefficients, 
m,  and  me  the  corresponding  moisture  concentrations  in  the  adjacent  air,  m,  ,ot  and 
m„t,ot  the  saturation  moisture  concentrations  corresponding  with  surface  temperatures 
T,  and  T„,  and  hv  is  the  latent  heat  of  vaporization  of  water.  According  to  Equation  22, 
Hm,  may  be  positive  or  negative,  as  determined  by  the  sign  of  the  bracketed  term.  Jt 
thus  represents  either  heating  by  condensation  or  cooling  by  evaporation.  Evidently, 
evaporative  cooling  can  only  occur  on  a  wet  or  icy  surface,  so  we  have  the  conditions 
on  Equation  22: 


Hm, ,  A  0  for  ws-1  =  0 
Hmt  e  A  0  for  w,  e  =  0 


where  w,_j  and  w,e  are  respectively  the  interior  and  exterior  surface  loadings  of 
water  or  ice. 

Hme|,  is  zero  except  during  melting  of  a  surface  ice  layer.  During  melting, 
the  surface  temperature  is  assumed  to  remain  constant  at  32°F.  The  appropriate 
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boundary  condition  is  then  chosen  from 

T,  =  32;  Tn  =  32  (24) 

to  replace  the  standard  boundary  condition  given  by  Equation  20.  The  requirement  for 
continuity  of  heat  flow  across  the  boundary  is  then  available  for  the  determination  of 
the  unknown  value  of  Hme,(.  Substituting  Equation  21  into  20  and  transposing  gives 

8T,  \ 

Hm.n  =  ki  9^7  +  h,(T,-32)  +  Hfh,  +  Hm,.l 

9Tn  {  (25) 

Hm.„  +  he(T„-32)  +  Hrh.e  +  Hmt.e 

for  the  inner  and  outer  surfaces  respectively. 


from: 


The  time  required  to  melt  a  given  surface  ice  loading  is  then  determined 

w,h,  =  f me"  Hm#l,dt  (26v 


where  w,  is  the  surface  ice  loading  and  h,  is  the  latent  heat  of  fusion  of  water.  The 
thermal  resistance  and  capacity  of  the  ice  layer  are  neglected  in  comparison  with 
those  of  the  windshield/backlight. 

The  surface  loadings  of  water  or  ice  are  determined  from  the  integrals: 

w,.,  =  /'  S,(mi_mi.s0t)dt  i  0  I 

(  (2?) 

wv.  *  jf*  Se(me-mn  l0,)dt  ^  0  J 

The  condition  w,  {  0  is  an  obvious  physical  limit. 

Illustration  of  Model  Characteristics 

The  general  characteristics  and  some  of  the  limitations  of  the  model  can  be 
more  fully  understood  by  reference  to  two  examples  in  which  the  behaviour  of  the  major 
variables,  namely  temperature  and  moisture  concentration,  is  followed.  These  same 
examples  are  also  used  in  later  sections  of  this  report  to  illustrate  the  use  of  the 
program. 

Example  1  -  Defrosting  of  a  Laminated  Windshield 


Figure  1  shows  plots  of  temperature  against  time  for  the  major 
points  of  interest  in  the  system. 

The  ambient  temperature,  and  initial  temperature  for  all  sub¬ 
systems,  is  20°F.  The  relative  humidity  is  100%  so  the  windshield  is 
assumed  to  be  coated  initially  with  0.01  oz/in2  of  ice,  other  effects  of  pre¬ 
cipitation  being  ignored.  The  external  airflow  over  the  vehicle  has  a  velocity 
of  one  ft/sec. 
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The  maximum  engine  coolant  temperature  is  limited  to  195°F  by 
the  thermostat.  Other  subsystem  parameters  are  also  assigned  typical 
values.  It  can  be  seen  in  the  figure  that  the  maximum  defroster  jet  temper¬ 
ature,  the  vehicle  interior  temperature  and  hence  the  local  defroster  jet 
temperature  are  directly  related  to  the  engine  coolant  temperature.  In 
practice,  while  the  coolant  temperature  rise  is  sharply  cut-off  by  thermostat 
opening,  some  further  rise  in  defroster  jet  temperature  may  occur.  The 
vehicle  interior  temperature  generally  continues  to  increase,  albeit  more 
slowly  than  at  first.  The  simplified  model  behaviour  results  mainly  from 
ignoring  the  thermal  capacities  of  the  heater  and  of  the  vehicle  interior. 

The  windshield  surface  temperatures  are  determined  primarily  by 
the  local  defroster  jet  temperature  which  is  a  function  of  distance  along  the 
jet  axis.  The  plotted  curves  refer  to  a  point  five  inches  from  the  base  of  the 
windshield  and  the  defroster  jet  slots.  The  effects  of  the  thermal  capacity 
and  resistance  of  the  laminated  windshield  which  are  represented  in  the 
model  can  be  seen  in  the  curves.  Also  of  interest  in  the  curve  for  the 
exterior  surface  are  the  small  discontinuity  at  32°F  associated  with  melting 
of  the  ice  layer  and  the  sharp  temperature  rise  following  the  completion  of 
evaporation  from  the  surface  at  28  minutes. 

In  the  present  example,  the  corresponding  curves  of  moisture 
concentration  are  of  no  special  interest.  Removal  of  exterior  ice  and  water 
by  sublimation  and  evaporation  occurs  throughout  the  first  28  minutes. 


Example  2  -  Befogging  of  an  Electrically  Heated  Backlight 

Figure  2  shows  curves  of  temperature  against  time  for  the  major 
points  of  interest  in  this  case. 

The  ambient  and  initial  temperature  is  33°F  with  100%  relative 
humidity.  The  backlight  is  therefore  cooled  by  a  film  of  rainwater.  The 
vehicle  has  six  occupants  with  wet  clothing.  The  engine  and  heater  charac¬ 
teristics  are  the  same  as  in  the  previous  example.  A  constant  heating  rate 
of  160  Btu/ft2hr  is  assumed  at  the  interior  surface  of  the  backlight,  which 
consists  of  a  single  5-inch  glass. 

Heat  transfer  between  the  vehicle  interior  and  the  backlight  occurs 
only  by  free  convection.  The  surface  temperature  rise  is  therefore  deter¬ 
mined  primarily  by  the  electrical  heating.  The  curves  refer  to  a  point  12.5 
inches  up  from  the  lower  edge  of  the  backlight  but  the  actual  defogging  per¬ 
formance  will  not  vary  markedly  with  position  in  this  case. 

Figure  3  shows  the  relevant  curves  of  moisture  concentration 
against  time.  The  mass  fraction  of  water  vapour  in  the  vehicle  interior  rises 
above  ambient  owing  to  the  vapour  generated  by  the  occupants’  breath  and 
clothing.  While  the  interior  vapour  mass  fraction  exceeds  the  saturation 
value  corresponding  with  the  backlight  interior  surface  temperature,  fogging 
occurs.  As  the  glass  warms  up,  the  surface  saturation  concentration  rises 
above  the  vehicle  interior  value  and  the  fog  begins  to  evaporate,  finally 
leaving  the  surface  dry  after  about  2|  minutes.  The  figure  illustrates  the 
rather  critical  dependence  of  the  predicted  defogging  performance  on  the 
assumed  characteristics  of  the  vehicle  interior  subsystem. 
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Discussion  of  the  Model 

An  ideal  model  of  an  automobile  defog/defrost  system  would  be  capable  of 
predicting  defog/defrost  contours  on  the  windshield/backlight  for  a  wide  range  of  system 
designs  and  operating  conditions.  The  present  model  clearly  falls  some  way  short  of 
this  ideal,  primarily  as  a  result  of  the  inadequate  level  of  detail  in  modelling  the  heat 
and  mass  transfer  processes  at  the  surfaces  of  the  windshield/backlight.  For  all  prac¬ 
tical  purposes,  the  requisite  level  of  detail  is  unattainable,  by  reason  of  the  very  large 
range  of  geometrical  and  other  parameters  which  determine  the  nature  of  the  flow  over 
the  surfaces.  It  is  further  apparent  that  so  long  as  the  distribution  of  local  heat  and 
mass  transfer  rates  can  not  be  predicted  with  any  confidence,  possible  refinements  in 
other  parts  of  the  model  may  offer  little  return  in  overall  accuracy  or  realism.  Thus, 
for  example,  consideration  of  in-plane  conduction  effects  in  the  windshield/backlight 
subsystem,  which  would  significantly  increase  complexity  and  computing  costs,  could 
not  be  justified. 


While  detailed  predictions  of  defog/defrost  contours  can  not  be  made  with 
the  existing  model  or  any  practicable  development  of  it,  its  value  is  nonetheless 
considerable.  It  may  be  expected  to  give  good  indications  of  the  effects  of  the  many 
available  system  parameters  on  the  defog/defrost  time  for  representative  locations  on 
the  windshield/backlight.  It  will  be  recalled  that  these  parameters  define  the  major 
characteristics  of  both  the  operating  conditions  and  the  vehicle  hardware.  Thus,  it  may 
be  used  for  preliminary  studies  of  alternative  performance  criteria  or  of  alternative 
system  designs. 

In  this  more  limited  context,  there  remain  several  parts  of  the  model  which 
could  usefully  be  improved.  The  improvements  rn.ght  include  the  following: 

(i)  extension  of  engine  and  heater  subsystem  model  to  include  engine  rpm, 

engine  load  and  coolant  mass  flow  effects; 

(ii)  substitution  of  a  differential  model  of  the  interior  heat  and  mass 

balances  for  the  existing  equilibrium  model; 

(iii)  development  of  a  representative  model  of  an  unheated  backlight 

defogger  blower; 

(iv)  the  facility  to  specify  different  start  times  for  certain  subsystems  or 

subsystem  components. 

From  a  preliminary  review  of  these  items,  it  appears  that:  (i)  could  be  represented 
by  additional  closed-form  expressions  derived  from  approximate  theory;  (ii)  would  in¬ 
volve  the  numerical  solution  of  a  simple  first  order  equation  and  would  slightly  increase 
computing  time;  (iii)  would  probably  require  some  experimental  studies  to  develop  a 
suitable  correlation  expression;  (iv)  seems  only  to  involve  program  modifications. 


USE  OF  THE  PROGRAM 
General  Description 


Figure  4  shows  the  major  features  of  the  program  in  flow  chart  form. 
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Each  input  data  set  comprises  one  set  of  vehicle  parameters  and  ambient 
conditions,  a  maximum  of  eight  locations  on  the  windshield/backlight  to  be  considered, 
and  a  specification  of  the  number  and  size  of  time  increments  to  be  used.  The  program 
has  a  simple  nested  loop  structure  in  which  the  time  history  of  the  ice  surface  tem¬ 
peratures  and  states  is  evaluated  for  each  specified  location  for  each  set  of  vehicle 
parameters  and  ambient  conditions.  The  surface  temperatures  and  states  are  written 
at  some  chosen  multiple  of  the  basic  time  increment. 


Specifying  the  Input  Data 

Table  I  prescribes  the  structure  of  the  input  data  set  which  is  oriented  to 
card  format  with  a  maximum  record  length  of  80  bytes.  The  first  field  starts  in  card 
column  1  {or  equivalent)  and  all  fields  are  contiguous  within  one  card  (record). 


The  following  notes  relate  to  Table  I: 

(1)  The  value  of  10 PT (!)  specifies  the  frequency  of  output  from  the  program. 
A  value  of  IOPT(l)  =  1  will  cause  surface  temperatures  and  states  to 
be  printed  at  every  time  increment,  a  value  IOPT(l)  =  2  at  every  other 
time  increment,  and  so  on. 

The  value  of  IOPT(2)  specifies  the  nature  of  the  interior  heat  transfer 
mechanism.  A  value  IOPT(2)  =  2  specifies  free  convection.  Any  other 
value  of  IOPT(2)  specifies  forced  convection  (the  defroster  jet). 

The  value  of  IOPT(3)  specifies  the  presence  or  absence  of  interlaminar 
heating.  A  value  IOPT(3)  =  0  causes  all  HRH(I)  to  be  taken  as  zero 
irrespective  c  the  values  read  in.  Any  other  value  causes  HRH(I)  to  be 
taken  as  read. 

10PT(4)  and  IOPT{5)  are  spare  and  should  be  left  blank. 

(2)  If  a  windshield  is  being  analyzed,  XO  is  the  distance  from  the  front  of  the 
hood  to  the  lower  edge  of  the  windshield.  If  a  backlight  is  being  analyzed, 
set  XO  equal  to  RFL. 

(3)  The  value  of  NDZ  must  be  specified  for  one  lamina  and  one  only.  This 
value  is  used  with  the  specified  time  increment  to  determine  a  'master* 
value  of  the  finite  difference  parameter  a  to  which  'slave'  values  of  a 
for  the  other  laminae  are  matched  as  closely  as  possible,  subject  to  other 
physical  and  program  constraints. 

(4)  When  IOPT(3)  =  0,  this  record  must  be  blank  (but  not  omitted)  for  cor¬ 
rect  program  functioning. 

(5)  The  next  record  may  be  either  a  data  set  delimiter  or  the  first  record  of 
a  new  data  set. 

Table  TT  shows  the  input  data  set  for  the  two  illustrative  examples  used  in  the 
present  report. 


Output  Format 


The  output  format  is  fixed  except  in  respect  of  the  choice  allowed  in  the 
value  of  IOPT(l).  Table  III  shows  the  program  output  for  the  two  examples. 

The  units  of  the  system  parameters  are  the  same  as  specified  in  Table  I  for 
input.  The  surface  state  integer  takes  the  values  one,  two  and  three  for  dry,  wet  and 
icy  surfaces  respectively. 

Run  Times  and  Stability 

Actual  run  times  in  any  particular  case  evidently  depend  on  a  number  of 
factors  some  of  which  are  under  the  control  of  the  program  user  while  others  are 
program-  and  installation-dependent. 

The  development  of  the  present  program  was  carried  out  on  an  IBM  360/67 
operating  under  TSS/360.  The  majority  of  examples  extended  over  180.0  increments  of 
10  seconds  (30  minutes  of  real  time).  Typical  CPU  times  for  these  conditions  were 
found  to  be: 

0.15  -  0.5  cpu  sec  to  read  data  and  set  up  problem; 

1.0  -  3.0  cpu  sec  to  compute  surface  temperatures  and  states  at  one 

location  over  a  30-minute  period. 

No  evidence  of  oscillation  or  divergence  due  to  rounding  errors  was  observed 
in  any  of  the  examples  used  in  program  development. 
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TABLE  I 


Preceding  page  blank 


ADSYM  INPUT  DATA  SET  STRUCTURE 


FTN  FIELD  FTN 

V'BLE  LENGTH  FORMAT 


NSTEP 


THETA 


TCMAX 

TCINF 


DESCRIPTION 


Problem  title,  preceded  by  0  (zero) 


Time  increment  in  finite  difference  solution 
Number  of  time  increments 
Number  of  locations  considered  (1  £  NX  £  8) 
Program  options111 


Initial  and  ambient  temperature 
Relative  humidity  percent 
Characteristic  length  of  car121 
Slant  height  of  windshield/backlight 
Wmdshield/baeklight  angle  to  horizon 
Length  of  car  roof 
Speed  of  car  (VCAR  >  0) 

Initial  coolant  temperature 
■Maximum  coolant  temperature 
Asymptotic  coolant  temperature  rise 
Characteristic  time  of  engine*  heater 
Heater  effectiveness  (0  s  I1EFF  £  1) 

Mass  flow  parameter  (0  £  WJT  s  1) 

Number  of  occupants 
Effective  mass  flow  through  car 

Width  of  defroster  jet  slot  (normal  to  windshield) 
Initial  jet  velocity 

Ratio  of  slot  length  to  width  of  windshield 

|  Number  of  laminae  s  NLAM  £  5) 
l 

|  Material  of  first  lamina 
|  Thickness  of  first  lamina 
Density  of  first  lamina 
Specific  heat  of  f.rst  lamina 
Thermal  conductivity  of  first  lamina 
Number  of  snblaminue13' 

Material  of  second  lamina 
Thickness  of  second  lamina 

etcetera 


Btu/lb”l 
Btu/hr  ' 


GO  I  6F1P.5  (NLAM*1)  values  of  interlaminar  heating  rates 


3J  8F10.5  (NX)  values  of  distance  along  reference  axis1 


TABLE  II 
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TABLE  III  (  cont'd  ) 
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APPENDIX  A 


RAINWATER  FILM  THICKNESS  ON  THE  WINDSHIELD/BACKLIGHT 


Basic  Equations 

Consider  a  rainwater  film  of  local  thickness  y  moving  with  local  mean  velo¬ 
city  v  down  the  inclined  surface  of  a  windshield/backlight,  as  shown  in  Figure  A-l. 
For  unit  width  of  the  film,  conservation  of  mass  within  the  elementary  length  ds  requires 
that: 

d(Pvy)  -  p  cos  9  dc  =  0  (A.  1) 

where  p  is  the  density  of  the  water,  p  the  precipitation  rate  and  9  the  inclination  of 
the  windshield/backlight  to  the  horizon. 

The  second  required  condition  for  the  determination  of  the  unknown  film 
thickness  and  mean  velocity  may  be  obtained  by  considering  the  equilibrium  of  the 
element  under  the  forces  acting  in  the  s  direction  due  to  its  inertia,  weight,  viscous 
drag  and  the  momentum  increment  due  to  precipitation.  Equilibrium  (or  conservation 
of  momentum)  requires  that: 

^  (Pvy)ds  -  gpy  sin  0  ds  +  r0ds  -  pvp  sin  9  cos  9  ds  «  0  (A.  2) 

where  t  is  time,  g  is  the  acceleration  due  to  gravity,  r0  is  the  viscous  shear  stress 
at  the  surface  and  vp  is  the  vertical  velocity  of  the  precipitation. 

Shear  Stress  t0 

The  shear  stress  r0  must  be  expressed  in  terms  of  the  local  film  thickness 
and  mean  velocity  before  (A.  2)  can  be  used.  It  is  assumed  that  the  film  flow  is  laminar 
and  characterised  by  a  simple  parabolic  velocity  profile.  Newton’s  expression  states 
that: 


where  p  is  the  dynamic  viscosity,  v*  is  the  velocity  at  distance  y'  from  the  surface. 
For  the  assumed  velocity  profile 

/dv*\  _  m0»  _  3v 

W’/y'  =  o  y  "  y 


3/iv 

y 


(A.  3) 


whence 


t. 
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Solution  of  Equations 


Substituting  for  r0  in  Equation  A.  2,  noting  that  ^  =  v  ^  and  that  p  is 


constant,  Equation  A.  2  becomes 


vp  ^  (vy)  =  gpy  sin  6  + 


pv  sin  6  cos  0  =  0 


But  from  Equation  A.  1  we  have 


so  that 


P  (vy)  =  P  cos  0 


vpcos0  -  gPy  sin®  +  -  pvp  sin0  cos0  =  0 


Consideration  of  the  relative  magnitudes  of  the  terms  in  this  equation  shows  that  the 
first  and  last  are  respectively  three  and  two  orders  smaller  than  the  remaining  terms. 
For  typical  windshield/backlight  angles  we  ■'an  therefore  write 

~  -  gpy  sin  0-0 

Multiplying  bv  y2,  rearranging  and  noting  v  =  p/p 


vy  - 


(A.  4) 


Equation  A.l  can  be  integrated  directly  to  give 


vy  =  ^  cos  9  +  (vy)4.0 


(A.  5) 


Equating  the  R.H.S.  of  Equation  A.  4  and  A.  5  we  obtain 


3  sin  0  sp  a  .  .  . 

HZ -  =  -f  cos  6  +  (vy)i=0 


whence 


{  3v 
}  sin 


cos  0 


(A.  6) 


The  value  of  (vy),.0  is  determined  by  assuming  that  one-quarter  of  the  pre¬ 
cipitation  falling  on  the  roof  flows  over  the  windshield/backlight,  that  is: 


(vy)i=0  * 


where  Lroo(  is  the  length  of  the  roof. 


PRECIPITATION 

RATE 

p  lb/ ft2  hr 


RAINWATER  FILM  MODEL 


-:«55??"'.'ir5 


•  ^  ,^^?*'v-v  rossa  - v > 
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APPENDIX  B 

FINITE  DIFFERENCE  SOLUTION  OF  WINDSHIELD/BACKLIGHT  SUBSYSTEM 

EQUATIONS 


Introduction 

Finite  difference  approximations  to  the  solution  of  partial  differential  equa¬ 
tions  are  sufficiently  familiar  that  no  great  elaboration  on  the  basic  approach  is  neces¬ 
sary  here.  Reference  7  provides  an  excellent  source  of  such  material.  The  following 
development  concentrates  on  the  special  features  of  the  present  problem. 

Basic  Equation  for  Interior  Temperatures 

Figure  B-l  shows  a  typical  grid  of  points  in  (z,  t)  co-ordinates  for  the  solu¬ 
tion  of  a  one-dimensional  heat  conduction  problem  in  a  laminated  material.  At  any 
interior  point  of  a  given  lamina,  the  heat  conduction  Equation  18  expressed  in  backward 
difference  form  is: 


(Az)2 


(B.  1) 


where  T^j  is  the  temperature  at  the  interior  point.  Letting  a  =  kAt./cP(Az)2  and 
rearranging  gives: 


T. 


oT,.,.,  +  (1  +  2a)Ti , 


aT 


i  +  5.) 


(B.2) 


where  it  will  be  noted  that  a  is  a  parameter  of  the  particular  lamina  considered. 


Continuity  Condition 


At  each  interface  between  laminae,  the  continuity  condition,  Equation  19 
may  be  used  to  determine  a  fictitious  value  of  T14.,  ,  denoted  T*1+, ,  for  use  in 
Equation  B.2.  In  finite  difference  form,  Equation  19  becomes: 


(TV 


i.j 


T-.,> 


Az, 


=  k- 


(Tl+1.,  -  T„) 
Az, 


whence,  T* +, , 


/  k2  Az,  \  k2  Az, 

V  "  k,  Az2)  T'-1  +  ki  Az2 


Az, 


T..,.,  +  (H,h),  X 


(B.3) 


where  the  significance  of  the  various  terms  is  apparent  from  Figure  B-l  and 
Equation  19.  It  will  be  noted  that  when  the  conductivity  and  grid  spacing  in  adjacent 
laminae  are  the  same  and  (Hth).  =  0,  the  required  result  for  a  homogeneous  material, 
namely  T*i+,  ,  =  T1+,  is  obtained. 


nil  —  ,-- 
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Substituting  Equation  B.  3  into  Equation  B.  2  gives  the  special  form  of  the 
basic  equation  for  points  on  an  interface: 


T.,-,  +  «,  —  (H,k), 


“.T-. ,  +  |!  +  (‘  *  £  3rj)j T 


k,  Az2/ 


^2 

a'  k,  Az2  r^,+,  i 


(B.  4) 


Boundary  Conditions 


The  boundary  conditions,  Equation  20,  may  similarly  be  used  to  determine 
fictitious  temperatures  T*_,  ,  and  T*14.,  ,  outside  the  boundaries  of  the  windshield/ 
backlight  for  use  in  Equation  B.  2.  In  finite  difference  form,  the  boundary  conditions 
become: 


(T„  -  T*.,.,) 


=  h.nt(T,,  -  T,nl)  -  (H,)im 


(T**,.,  -  T.  ,) 

K  - AS -  =  "  Tm>  +  <Hs)„, 


The  significance  of  the  various  terms  should  again  be  clear  from  Equation  20  and 
previous  development  in  this  appendix. 

Transposing  the  above  equations,  we  obtain 

/  h,n,AzA  Az,  l  -I 

=  v 1  -  -kH  t-,  ^  —  lh-T- +  <h.>J 

;  (B.5) 

/  h«„Azn\  Azn  ,  )\ 

T*..u  =  -  -X~)  T,.,  +  Jh„,Text  +  (Hi)ext; 

Substituting  Equation  B.5  into  the  basic  equation  gives  the  special  forms  associated 
respectively  with  the  interior  and  exterior  surfaces  of  the  windshield/backlight: 


+  tt  ih-T-  +  <»>>„!  -  1  +  <*.  f1  +  h-’  — n T'.. +  a-T-- 


T,,  +  a.T,.,. 


(B.6) 


+ «-  —  +  +  1  +  m1  +  h-  tt, 
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Initial  Condition 

The  basic  Equation  B.  2  for  interior  points  w:‘h  the  special  forms,  Equa¬ 
tions  B.4  and  B.6,  for  points  on  interfaces  and  boundaries  provides  a  sufficient  number 
of  equations  for  the  determination  of  the  current  temperature  distribution  T,  from 
the  previous  temperature  distribution,  T,  , . 

The  initial  condition: 

T,  ,  =  T0  for  all  i  (B.7) 

completes  the  data  required  to  determine  the  temperature  distribution  for  any  j  >  1. 
Ta  is  a  specified  uniform  initial  temperature,  equal  to  the  ambient  temperature  in  the 
present  case. 

Melting  of  Surface  Ice 


Equation  B.6  for  the  surface  temperatures  is  invalid  during  melting  of 
surface  ice  films.  However,  owing  to  the  backward  difference  formulation  of  the 
equations,  it  will  not  be  apparent  that  one  or  other  has  become  invalid  until  after  it  has 
been  used  to  determine  a  surface  temperature  which  exceeds  the  melting  point  cf  ice. 
To  avoid  the  misapplication  of  Equation  B.6  and  the  consequent  necessity  to  recalculate 
temperature  distributions  following  the  occurrence  of  melting,  an  estimate  of  the  final 
surface  temperatures  is  required  before  Equation  B.6  is  applied. 

Such  an  estimate  is  available  through  the  use  of  the  forward  difference 
approximation  to  the  governing  equations.  In  forward  difference  form,  we  have  instead 
of  Equation  B.  1: 


(B.8) 


and  instead  of  Equation  B.2: 


Mi  =  aT,,.,  +  (1  -  2a)T,  ,  +  o-T,^, 


(B.  9) 


The  boundary  conditions  expressed  in  Equation  B.5  are  unchanged,  so  sub¬ 
stituting  these  in  turn  into  Equation  B.9  leads  to  the  following  estimates  of  the  surface 
temperatures  based  on  the  forward  difference  approximation: 


^zi  <  •  i  (  /  ^zi' 

T„.,  -  «  —  *  (H,)mt;  *  jl  -  or  ^1  ♦  h.„,  — 


rr  ,  ) 


-  «  —  Jh...T.„  +  (Hm)extj  Ml  -  a  (1  +  h„, 


T  +  aT  , 


Tm  +  aT,-.., 


(B.  10) 
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In  the  event  that  one  or  other  of  these  equations  predicts  that  the  temperature 
of  an  ice-covered  surface  will  exceed  32°F  during  the  next  time  increment,  the  alter¬ 
native  boundary  condition: 


=  32 


(B.  11) 


is  applied  to  either  or  both  of  the  first  and  last  equations,  reducing  the  number  of 
unknown  T,  (  accordingly. 

The  finite  difference  form  of  Equation  25,  namely 


H 


melt.mt 


H 


melt. ext 


kn 

Az“  (T-U 


-  T,,)  +  hint(Tint 

-  T(  |)  +  h„,(T„, 


-  T;  , )  +  (H,h)t  +  (Hmt)\ 

l  (B.  12) 

-  T.  ,  )  +  (H,h)(  +  <Hmt)J 


provides  the  basis  for  determining  the  heat  available  for  melting  the  surface  ice  layer 
during  the  j'h  time  increment. 

The  time  required  to  melt  the  surface  ice  layer  is  obtained  from  the  finite 
difference  form  of  Equation  26: 


(Wjh()  =  £  H  melt  mell 


(B.  13) 


The  surface  loading  ws  is  determined  from  the  finite  difference  form  of 
Equation  27: 


=  Z  g.M,  ~  m,  )  At  0 
1 
J 

=  Z  Sex..,  At  {  0 


(B.  14) 


which  apply  quite  generally  to  ice  or  water  and  irrespective  of  the  occurrence  of 
melting. 

In  general,  Atmel,  in  Equation  B.  13  is  simply  At.  However  in  the  first  time 
increment  in  which  melting  occurs,  Atme„  is  determined  from: 


Atmelt  T,  |+,  -  32 
At  =  Tm*.  -  T„ 


where  T,  is  the  forward  difference  estimate  of 
the  surface  temperature  in  the  absence  of  melting. 
In  the  last  increment  in  which  melting  occurs, 
Atmel)  is  chosen  to  be  just  large  enough  to  satisfy 
Equation  B.  13. 


A  ~ J  •■*:*:  **  ‘■'i  *  .  A*  . 


T'K-S'AV'-soi^  -  i;j. 


INTERIOR 

SURFACE 


INTERFACES 


EXTERIOR 

/SURFACE 


BHi 


H 


SCHEMATIC 

DISTRIBUTION 


FIG.  B-l  :  TYPICAL  GRID  FOR  FINITE  DIFFERENCE  APPROXIMATION 
T0  WINDSHIELD/ BACKLIGHT  TEMPERATURE  DISTRIBUTION 


kVV3WRy** y—,*^ f,x  ,  ;^t,'  .,' ,7l^^5>“;vvr. 


Preceding  page  Plan 


C  NRC/NAE  AUTOMOBILE  DEFOG/DEFROST  SYSTEM  STUDY 
C  HEAT  AND  MASS  TRANSFER  MODEL 

C  MAIN  PROGRAM  VERSION  1  LEVEL  2  DECEMBER  1971  ERW 
C  MODULE  NAHE=ADSYM 

DIMENSION  AMAT  (2, 5)  ,  X  (8)  ,IOPT{5) 

COMMON/ONE/TO , P HUM, XO, SLHT, THETA ,  P.FL,  VCAR ,  TCO,TCM AX, TCINF,  TAU, 
&HEFF, WJT,WJM , ZJET, VJET, SPREAD  ,  NOCC/TKO/THK  (5)  , DENS  (5)  ,  SPHT  (5)  , 
5COND  (5)  ,  ALPHA  (5)  ,DZ{5)  ,H?.H  (6)  ,IIDZ  (5)  , NLAM,  NIF,NGP/’rHP.EE/A  (3,51)  , 
STOLD (51)  ,TNEW  (51) 

2  READ  (1,9 99, EM D= 100) 

READ  (1,997)  DTIM,NSTEP,  NX,  (IOPT(I)  ,1=1,5) 

READ  (1,995) TO , RHUM, XO, SLUT, THETA ,RFL, VCAR 

READ  (1,995) TCO,TCMAX,TCINF,TAU, HEFF, WJT 

READ ( 1,^92) NOCC, WJM 

READ  (1,995) ZJET, VJET , SPREAD 

READ (1,990)  NLAM 

READ  (1,987)  ((AMAT(J,I)  ,J  =  1,2)  ,THK  (I)  ,DENS(I)  ,SPHT  (I)  ,COND(I)  , 

&NDZ (I) ,1=1, NLAM) 

NIF=NLAM+1 

RE  AD  ( 1 ,995)  (HRH  (I)  ,1=1  ,NIF) 

READ (1,995)  (X  (I) ,1=1 , NX) 

CALL  STEPS (DTIM) 

WRITE  (3,985) 

WRITE  (3,999) 

WRITE  (3,982) TO, RHUM, XO ,SLHT, THETA, RFL, VCAR 
WRITE  (3,980) TCO,TCMAX,TCINF, TAU, HEFF, WJT 
WRITE  (3,977)  NOCC,WJM 
WRITE  (3,975) ZJET, VJET, SPREAD 
WRITE  (3,972) 

DO  5  1=1, NLAM 
WRIT  3  (3 , 970)  HRH  (I) 

THKP=12. *THK (1) 

5  WRITE  (3 , 96V)  j,  (AMAT(J,I)  ,J-1,2)  ,THKP,NDZ(I)  , DENS  (I)  ,  SPHT  (I)  , 

&COND  (I)  ,  ALPHA  (I) 

WRITE  (3, 97C  HRH  (NIF) 

CALL  DEWPNT vTO, RHUM, TD2 , VMF2) 

CALL  SETUP 
DO  15  1=1, NX 
WRITE  (3,965)  X  (I) 

CALL  EXTSS (I, X  (I) ,TF2, HTC2,TCM2) 

TIM=D 

ISTE.  0 

KPRNT=IOPT (1)  -1 
DO  7  J=1 ,NGP 
7  TNEW  (J)  =TO 

10  CALL  EHSS (TIM ,TC,TJET) 

CALL  INTSS (TIM,TJET, VMF2 ,TCAR,VMFCAR) 

CA'L  I3LSS (IOPT  (2) ,I,X(I)  ,TIM,TJET,VKF2,TCAR,VMrCAR,T?1, VMF1 ,HTC1, 


CALL 

CA'L 

STCM1) 

CALL 

CALL 


CALL  INSCON  (TIM, DTIM, TF1  ,VMF1  ,i.  ;c  1  ,TCM  1 ,  KS 1 ,  WS1 ,  H SI ,  Ml) 
CALL  EXSCON (TIM, DTIM, TF2 , VMF2 ,HTC2,TCM2,KS 2, WS2,HS2, M2) 
KPRNT=KPRNT+ 1 

IF(KPRNT.NE.IOFT(1))  GO  TO  12 

WRITE  (3,962)  TIM, TNEW  (1)  ,KSl,TNEW  (NGP)  ,  KS2 


ts D t,  inti, ,  y. «  _c 


i, 

F 


-v 

fv 

t 

if 


r 


\ 

i' 

f 

;■ 

j 

« 


I' 


KPRNT- 0 

12  ISTEP=ISTEP+ 1 

IF (ISTEP.GT. NSTEF)  GO  TO  15 
TIM=TIM+DTIM 

CALL  MODEQ (IOPT  (3) ,T?1 ,HTC 1 , HS 1 , M ,TF2 ,HTC2, HS2, K 2, 1 1, 12) 

CALL  GE06 (II, 12, A, TN EH, TOLD) 

GO  TO  10 
15  CONTINUE 
GO  TO  2 
100  STOP 

999  FORMAT  ('JOB  TITLE  NOT  EXCEEDING  56  CHARACTERS  ») 

997  FORMAT(F10.5,7I5) 

995  FORMAT (8F1 0. 5) 

992  FORMAT (15, F10. 5) 

990  FORMAT  (15) 

987  FORM AT (2A4,4?10. 5,15) 

985  FORMAT{« 1AUTOMOBILE  DEFOG/DEFROST  SYSTEM  MODEL  (ADSYM)  •) 

982  FORMAT (‘OEXTERIOR  SUBSYSTEM  PARAMETERS'/  5X, 'AMBIENT  TEMI ERATUFE  = 
S' ,F6. 1,SX, 'RELATIVE  HUMIDITY  =',F6.1,5X, '  'HA 3ACTE FISTIC  LENGTH  OF 
SCAR  =',F5.  1/5X, 'SLANT  HEIGHT  OF  WINDSKIELD/BACI'7  IGHT  = ' ,  F5.  1 , 5X, '  S 
SLOPE  (DEG)  =• ,F5.1,5X, 'LENGTH  OF  ROOF  =• , F5. 2,5X , • SPEED  =',F6.1) 
980  FORMAT ('OENGINE  AND  HEATER  SUBSYSTEM  PARAMETERS' /5X, • INITIAL  COOLA 
SNT  TEMPERATURE  =', F6 . 1 , 5X, • MAXIMUM  COOLANT  TEMPEFATURE  =',F6.1,5X, 
S  '  ASYMPTOTIC  COOLANT  TEMPERATURE  =  *  , F 6 .  1/5X ,  •  CHAR ACTEP.I STIC  TIME  =' 
S,F€. 1,5X,' HEATER  EFFECTIVENESS  = ' , ?5. 2 , 5X, • MASS  FLOW  PARAMETER  =• , 
ST5.2) 

977  FORMAT  (' OINTEF.IOP.  SUBSYSTEM  ?ARAMETEES'/5X,  •  NUMBER  OF  OCCUPANTS  =• 
S, 13, 5X, 'EFFECTIVE  MASS  FLOW  =',F7.1) 

975  FORMAT('ODEFRCSTER  JET  SUBSYSTEM  PARAM ET ERS ' /5X , • SLOT  WIDTH  =',F5. 

S2,5X, 'INITIAL  JET  VELOCITY  =',F6. 1,5X, 'SPREAD  =',F6.3) 

972  FORMAT(»  OWINDSHIELD/BACKLIGHT  SUBSYSTEM  PA RA METEF S' /8X ,' LAMINA ’, 6X 
S' MATERIAL*  ,5X,  '  IHICKNESS ' , 5X , '  NUMBER  OF' ,  7X ,' DFNS ITY ',  6X  ,'  SPECIFIC¬ 
S'  , 2X, 'CONDUCTIVITY ', 9X, » ALPHA », NX, ' HEATING */46X , ' SUBLAMINAE' , 24X, • 
SHEAT' ,38X, 'RATE') 

970  FORMAT(112X,F14.2) 

967  FORMAT(I14,6X,2A4  ,F14. 3,1 14,?14 . 2, 2F1 4. 3 ,F1 4. 2) 

965  FORMAT ( • 1PREDICTED  SURFACE  CONDITIONS  AT»,F5.1,»  INCHES  FROM  BASE 
SOF  WINDSHIELD/EACKLIGHT'/26X, 'INTERIOR  SURFACE' , 1 2X, • EXTERIOR  SURF 
S ACE* /I OX,' TIME ' ,8X,' TEMPERATURE • ,4X, ' STATE' , 8X , 'TEMPERATURE* ,4X ,' S 
STATE') 

962  FORMAT (8X,F6. 1, 14X,F5. 1 , 8X ,1 1 , 14X,F5. 1 ,8X ,11) 

END 
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SUBROUTINE  STEPS (DTIM) 

C  DETERMINES  ALPHA  FOR  FINITE  DIFFERENCE  EQUATIONS 
C  AND  COMPUTES  NUMBER  OF  GRID  INTERVALS 
C  MODULE  NAME=  ADSM1 

COMMON/THO/THK  (5)  #  DENS  (5)  ,  SPHT(5)  ,COND  (5)  ,  ALPHA  (5)  ,DZ(5)  ,HRH(6)  , 
SNDZ  (5)  , NLAM, NIF , NGP 
DO  2  I=1,NLAM 
THK(I) =  . 08333  333*THK (I) 

IF (NDZ (I) « NE. 0)  IS=I 
2  CONTINUE 

DT=DTIM*. 2777778E-03 

DZ  (IS)  =THK  (IS)  /FLOAT  (NDZ  (IS) ) 

ALPHA  (IS)  =COND  (IS)  *DT/ (SPHT  (IS)  *DENS  (IS)  *DZ(IS)  *DZ(IS)  ) 

NGP=  1 

DO  7  I=1,NLAM 
IF (I. EQ. IS)  GO  10  5 
CONST=SPHT (I) *DENS  (I) 

CONST=COND  (I) /CONST 

DZ  (I)  =SQP.T  (CONST*DT/ALPH  A  (IS)  ) 

NDZ (I)  =INT  (THE  (I)/DZ  (I) ) 

IF  (NDZ (I) .LT. 2)  NDZ (I) =2 
IP(NDZ(I)  .GT.  10)  NDZ  (I)  =10 
DZ  (I)  =THK  (I)  /FLOAT  (NDZ  (I) ) 

ALPHA (I) =CONST*DT/ (DZ (I) *DZ (I) ) 

5  NGP=NGP*NDZ(I) 

7  CONTINUE 
RETURN 
END 


SUBROUTINE  DEWFNT  (T.RH,TD,VMF) 

C  COMPUTES  DEW  POINT  AND  MASS  FRACTION  OF  WATER  VAPOUR  IN  MOIST  AIR 
C  FROM  THE  MAGNUS  FORMULA.  REFERENCE:  IHVE  HANDBOOK 
C  MODULE  NAME=ADSM 2 

TC=. 5555555* (T-32.) 

IF  (T.LT.  32.)  GO  TO  2 
PLG=7.5*TC/(237.6+TC) +.78571 
GO  TO  5 

2  PLG=9. 5*TC/( 265. 5+TC) +.73571 
5  PVS= 1 0 . **PLG 
PV=RH*PVS*.0 1 
PLG=ALOG 10  (PV)-. 78571 
TD=237. 3*PLG/  (7. 5-PLG) 

IF (TD.LT.O.)  TD=265.5*PLG/(9. 5-PLG) 

TD=TD* 1.8+32. 

PAIP.=  10 1 3.  25-  PV 
VMF=18.*PV 

VMF=VMF/ (29.+PAIR+VMF) 

RETURN 

END 


s 
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SUBROUTINE  SETUP 

C  FORMS  MATRIX  OF  COEFFICIENTS  IN  FINITE  DIFFERENCE  EQUATIONS 
C  FOR  WINDSHIELD/BACKLIGHT  TEMPERATURE  DISTRIBUTIONS 
C  MODULE  NAME=ADSM3 

COMMON/TWO/TH K  (5)  ,DENS  (5)  ,  SPHT(5)  ,COND  (5)  ,  ALPHA  (5)  ,DZ(5)  ,HRH(6)  , 
SNDZ(5> ,NLAH, NIF,NGP/THREE/A(3,51) ,TOLD(51)  ,TNEW  (51) 

A(1,1)=0. 

A(3,  1)=- ALPHA  (1) 

12=1 

DO  5  IL= 1, NLAM 
11=12*1 

12=1 1+NDZ (IL) -2 
FA=1 . *2. *ALPH A  (IL) 

DO  2  1=11,12 
A  (1,I)=-ALPHA  (IL) 

A  (2, 1)  =FA 
2  A  (3,1)  =A  (1 ,1) 

12=12+1 

IF  (IL.GE. NLAM)  GO  TO  7 

DZK=DZ  (IL+ 1)  *COND  (IL)  /  (DZ  (IL)  *COND  (IL+ 1) ) 

A  ( 1 , 12)  =-ALPH  A  (IL) 

A  (2, 12)  =  1 .  +ALPHA  (IL)  *  ( 1 .  +DZK) 

5  A  (3,12)  =-ALPHA  (IL)  *DZK 
7  A  (1,12)  =-ALPH  A  (NLAM) 

A  (3,12) 

RETURN 

END 


is 
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SUBROUTINE  EXTSS (I,X,TF,HTC,TCM) 

C  EXTERIOR  SUBSYSTEM  MODEL 

C  COMPUTES  HEAT  AND  MASS  TRANSFER  COEFFICIENTS  ON  EXTERNAL  SURFACE  OF 
C  WINDSHIELD/BACKLIGHT  UNDER  SPECIFIED  AMBIENT  CONDITIONS 
C  MODULE  NAME=ADSM4 

COMMON/ONE/TO, RHUM, XO, SLHT, THETA, RFL, VCAR,TCO,TCMAX,TCINF,TAU, 
SHEFF,WJT,WJM,ZJET,VJET,SPREAD,NOCC 
IF  (I.GT.  1)  GO  TO  5 
TF=TO 

IF (TO. GE. 25. . AND. RHUM. GE. 100.)  GO  TO  2 
C**ESTABLISH  CONSTANTS  FOR  ZERO  PRECIPITATION  OR  DRY  SNOW** 

GNU=. 135E-03 
COND=. 135E-01 
PR=. 710 

FPR=EXP (.3 3* A LOG  (PR) ) 

GO  TO  7 

C**ESTABLISH  CONSTANTS  FOR  HEAVY  RAINFALL** 

2  PI=3. 141593 

TH=THETA*PI/1 80. 

CTH=COS(TH) 

PDOT=5. 205 
GNU=. 6930E-0 1 
RHO=62. 57 
G=32. 18*3600. 

COND=- 345 

CONST=3. *GNU*PDOT/ (G*RHO*SIN {TH) ) 

GO  TO  10 

5  IF (TO.GT.25. .AND.RHUM.GE.100.)  GO  TO  10 
C**COMPUTE  LOCAL  HEAT  AND  MASS  TRANSFER  COEFFICIENTS  FOR  DRY  CONDITIONS* 
7  XL=XO*. 083333  33 *X 
RE=VCAR*XL/GNU 
FRE=EXP(.8*ALOG{RE)) 

HTC=.0296*FPR*FRE/XL 

TCM=4.8*HTC 

RETURN 

C**COMPUTE  LOCAL  HEAT  TRANSFER  COEFFICIENT  IN  HEAVY  RAINFALL** 

10  XL=. 08333333*  (SLHT-X) 

H=EXP (. 3333333*ALOG (CONST* (XL*CTH  +  . 25*RFL)  ) ) 

HTC=COND/H 

TCM=0. 

RETURN 

END 
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SUBROUTINE  EHSS  (TIN#TC#TJET) 

C  ENGINE  AND  HEATER  SUBSYSTEM  MODEL 

C  COMPUTES  COOLANT  AND  DEFROSTER  JET  TEMPERATURES  AT  SPECIFIED  TIME 
C  MODULE  NAME=ADSM5 

COMMON/ONE/TO, RHUM,XO,SLHT, THETA , RFL, VCAR,TCO#TCM AX, TCINF,TAU, 
SHEFF,WJT,WJM# Z JET, VJET, SPREAD, NOCC 
IF (TIM.GT. 0. )  GO  TO  2 
TCLAST=0. 

2  IF (TCLAST. GE. TCMAX)  GO  TO  5 

TC=TCO+ (TCINF-TCC) * ( 1 . -EXP (-TIM/T AU)  ) 

IF  (TC.GT. TCMAX)  TC=TCMAX 
TJET=HBFF* (TC-TO) +TO 
TCIAST=TC 
TJLAST=TJ£T 
RETURN 
5  TC=TCLAST 
TJET-TJLAST 
RETURN 
END 
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SUBROUTINE  INTSS(TIM#TJET,VMF2,TCAR,VMFCAR) 

C  VEHICLE  INTERIOR  SUBSYSTEM  MODEL 

C  COMPUTES  INTERIOR  AIR  TEMPERATURE  AND  MOISTURE  CONTFNT  AT  SPECIFIED  T 
C  MODULE  NAM£=ADSM6 

COMMON/ONE/TO, RHUM, XO, SLHT, THETA , RFL, VCAR,TCO,TCMAX,TCINF,TAU, 
GHEFF, WJT, WJM,Z JET, VJET, SPREAD, NOCC 
IF  (TIM. GT. 0. )  GO  TO  2 
GOCC=10. 

IF  (RHUM. LT. 100.)  GOCC  =  0. 

ANO= FLOAT (NOCC) 

2  TCAR=TO+WJT* (TJET-TO) 

CALL  DEWPNT (TCAR, 100. , TD,VMFSAT) 

VMFCAR=VMF2*ANO* (. 16+GOCC* (VMFSAT-VMF2) ) /(WJM+ANO*GOCC) 

RETURN 

END 


SUBROUTINE  IBLSS  <IOPT,I,X,TIM,TJET,  VMF2,TCAR,VMFCAR,TP'*,  VMF1.HTC, 
STCM) 

C  INTERIOR  BOUNDARY  LAYER  SUBSYSTEM  MODEL 

C  COMPUTES  TEMPERATURE,  MOISTURE  CONTENT,  LOCAL  HEAT  AND  MASS 
C  TRANSFER  COEFFICIENTS  FOR  INTERIOR  BOUNDARY  LAYER 
C  MODULE  NAME=ADSM7 

COMMON/ONE/TO, RHUM, XO, SLHT, THETA, RFL, VCAR, TCO.TCM AX, TCINF,TAU, 
8HEFF,WJT, WJM,ZJET,VJET, SPREAD, NOCC/THREE/A (3,51) ,TOLD(51)  ,TNEW(51) 
IF  (IOPT. EQ. 2)  GO  TO  7 

C**ESTABLISH  CONSTANTS  FOR  FORCED  CONVECTION** 

IF(I.GT.I)  GO  TC  2 
A  1=.  5222E-06 
A2=~. 4938E-09 
A 2=-. 7 1604 E- 09 
Bl=. 2544E-04 
B2=— . 6172E-08 

ZEFF=ZJET*SPREAD*. 08333333 
FZE?F=EXP (. 325*ALOG (ZEFF) ) 

2  IF  (TIM.GT. 0. )  GO  TO  5 
XFT=. 08333333*X 
XPFM=XFT*6.666667*ZEFF 
FXPRM=EXP (-ALOG (XPRM)  ) 

EECF=SQRT ( 1. 6/  <  1 .  ♦ .  1 5*XFT/ZEFF) ) 

IF  (DECF « GT. 1 . )  DECF=1. 

C**COMPUTE  HEAT  AND  HASS  TRANSFER  COEFFICIENTS  FOR  FORCED  CONVU.  rxON** 

5  FTJ1=TJET*10. 

FTJ2=FTJ1* (TJET-80.) 

GNU=. 1288E-03+A1*FTJ1+A2*FTJ2 
COND=. 1287E-0 1 ♦B1*FTJ1*B2*FTJ2 
TF 1=TCAR*DECF*  (TJET-TCAR) 

VMF1=VMFCAR+DECF* (VMF2-VMFCAR) 

HTC=.  36*COND*EXP  (.65*ALOG  (VJET/GNU) )  *FZEFF*FXPF.M 

TCM=4.8*HTC 

RETURN 

C**COMPUTE  HEAT  AND  MASS  TRANSFER  COEFFICIENTS  FOR  FREE  CONVECTION** 

7  TF1=TCAR 
VMF1=VNFCAR 

TDIFF=ABS (TCAR-TNEH (1)  ) 


IF  (TDIFF) 10,  10,  12 


10  HTC=0. 

GO  TO  15 

12  HTC=.205*EXP(.33*ALOG (TDIFF)) 
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SUBROUTINE  INSCON (TIM , DTIM ,TF, VMF, HTC, TCM, KS , WS, HS. MELT) 

C  DETERMINES  INTERIOR  SURFACE  CONDITION  OF  WINDSHIELD/BACKLIGHT 
C  AND  COMPUTES  TOTAL  SURFACE  HEATING  RATE  FOR  NEXT  STEP 
C  MODULE  NAME=ADSM8 

COMM ON/ONE /TO, RHUM, XO, SLHT , THETA , RFL , VCAR, TCO, TCM AX, TCINF, TAU, 
SHBFF, WJT, W JM , ZJET, VJFT,SPREAD , NOCC/TWO/THK (5) , DENS (5) , SPHT (5)  , 
SCOND  (5)  ,  ALPHA  (5)  ,DZ(5)  ,  HRH  (6)  ,  NDZ  (5)  ,  NLAM,  NIF,  NGP/THREE/A  (3,51)  , 
&TOLD  (51)  ,TNEW  (51) 

IF (TIM. GT. 0. )  GO  TO  5 
DT=DTIM*.2777778E-03 
CDZ=COND  (1)  /DZ  (1) 

ADC=ALPHA ( 1) /CDZ 

C**INITIALIZE  SURFACE  CONDITION** 

MELT=0 
KS=  1 
WS=0. 

GO  TO  10 

C**DETEPMINE  CURRENT  SURFACE  CONDITION** 

5  IF (WS . GT  .0 .)  GO  TO  7 
KS=  1 

GO  TO  10 
7  KS=2 

IF  (TNEW  ( 1)  .LE.  32.)  KS=3 

C**COK PUTE  TOTAL  HEAT  AND  MASS  TRANSFER  DURING  NEXT  STEP*- 
10  CALL  DEWPNT (TNEW  (1) , 100. ,TDP,VMFSAT) 

DWS=TCM*  (VMF-VMFSAT) 

IF  (DWS)  12,  15,  15 
12  IF  (WS . LE. 0 . )  DWS=0. 

15  WS=WS+DWS*DT 

IF(TNEW(1) -32.)  17,  17,20 
17  HS=HP.H  (1)  +DWS*1204. 1 
GO  TO  22 

20  HS=HRH(1) +DWS*1060.61 
22  IF  (KS . NE. 3)  RETURN 

c**TFST  FOR  OCCURRENCE  OF  MELTING  DURING  NEXT  STEP** 

IF  (MELT.GT.O)  GO  TO  25 

TNEXT=ADC* (HTC*T?+HS) + (1 .-ALPHA (1)-HTC*ADC) *TNEK ( 1) + ALPHA ( 1) * 
6TNEW  (2) 

IF (TNEXT.LT. 32. )  RETURN 

C**COM?UTE  MELTING  RATE  AND  TEST  FOR  COMPLETION** 

MELT= 1 
HMSUM=0. 

DTM=DT* (TNEXT-32.)/(TNEXT-TNEW (1)  ) 

GO  TO  27 
25  DTM=DT 

27  HMSUM=HMSUM+DTM* (CDZ* (TNEW (2)  -TNEW (1) ) +HTC* (TF-TNEW ( 1 )  )+HS) 
HMREQ=143.49*WS 
DEF=HMREQ-HMSUM 
IF  (DEF)  30, 30,32 
30  MELT=0 

HS=- DEF/DT 
32  RETURN 
END 
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SUBPOUTIME  EX  SCON (TI M, DTIM ,TF , VKF, HTC, TCM , KS , WS , H S, MELT) 

C  DETERMINES  EXTERIOR  SURFACE  CONDITION  OF  WI ND SH IELD/B ACKL IG HT 
C  AND  COMPUTES  TOTAL  SURFACE  HEATING  RATE  FOR  NEXT  STEP 
C  MODULE  NAME=ADSM9 

COMKON/ONS/TO,RHUM,XO,SLHT, THETA,  RFL  ,  VCAP. ,  TCO,  TCP  AX,  TC  IN  F  ,  TAU  , 

?■  H EFF,  W JT, W  JM,  2 J 2T,  V JET,  SPREAD,  NOCC/TWO/TiiK  (5)  ,D?HS  (5)  ,  SPHT  (5)  , 
8COND  (5)  ,  ALPHA  (5)  ,  DZ  (5)  ,HP.H{6)  ,ND2  (5)  ,  NLA  M  ,  N1F,  NGP/THEF,  E/A  (  3 , 5  1 )  , 
&TOLD  (51)  ,TNEW  (51) 

IF  (TIM.GT.  0.)  GO  TO  5 
DT=DTIM*.2777778E-03 
CDZ=COND (NLA  M) /DZ (NLAM) 

ADC= ALPH  A (NLA  M) /CDZ 
C**INITIALIZE  SURFACE  CONDITION** 

M  E  LT  =  0 

IF  (R  HUM . GE . 100.)  GO  TO  2 
KS=  1 
WS  =  0. 

GO  TO  10 
2  KS=2 
WS=. 09 

I F  (TO . LT .  ) 2. )  i\S=3 
GO  TO  10 

C**DETERMINS  CURRENT  SURFACE  CONDITION** 

5  IF  (WS.GT.O.)  GO  TO  7 
KS=  1 

GO  TO  10 
1  KS=2 

IF(TNEW  (NGP)  .  LE.  32.)  KS=3 

C**COM PUTE  TOTAL  HEAT  AND  MASS  TRANSFFR  DURING  NEXT  STEP** 

10  CALL  DEVPNT (TNEW  (NG?)  , 100. ,TDP, VMFSAT) 

DWS=TCM*  (VMF-  VMFSAT) 

IF  (DNS)  12,  15,  15 
12  IF  (WS.IE.O.)  DWS=0. 

15  WS  =  WS*DWS*  DT 

IF  (TNEW  (NG?)  -32.)  17,  17,20 
17  HS  =  HRH(NIF) +DWS*1204.  1 
GO  TO  22 

20  HS=HRH(NIF) +  DWS* 1 060 . 6  1 
22  I?(KS.NF.3)  RETURN 

C**TFST  FOR  OCCURRENCE  OF  MELTING  DURING  NEXT  STEP** 

IF  (MELT. GT . 0)  GO  TO  25 

TNEXT=ADC*  (HTC*TF+HS)  -  ( 1 . -ALPHA  ( NLAM)  - HTC* ADC)  *TNEW  (NGP)  + 

S ALPHA  (NLAM)  *TNEW  (NGP-  1 ) 

IF  (TNEXT.LT. 32.)  RETURN 

C**COMPUTE  MELTING  RATE  AND  TEST  FOR  COMPLETION** 

MELT=  1 
HMSUM=0. 

DTM=DT* (TNEXT-32.)/(TN EXT- TNEW (NGP) ) 

GO  TO  27 
25  DTM=DT 

27  HMSUM  =  HMSU  MOTM*  (CDZ*  (TNEW  (NGP-  1 )  -"’NEW  (NGP)  )  *HTC*  (TF-TNEW  (NGP)  )  + 
SHS) 

HMREQ=143. 49*WS 
DEF=HMREC~HMSUM 
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SUBROUTINE  MOD 2Q <K,TF1 , HTC 1, HS1 , M1,TF2 ,HTC2, HS2, M2, II , 12) 

C  MODIFIES  FINITE  DIFFERENCE  EQUATIONS  TO  INCLUDE  CURRENT  VALUES 
C  OF  TIME  DEPENDENT  QUANTITIES 
C  MODULE  NAME=ADSM1 0 

COMMON/TWO/THK (5)  ,  DENS  (5)  ,SPHT(5)  #COND  (5)  #  ALPHA  (5)  ,DZ(5)  ,HRH  (6)  , 
&NDZ  (5) , NLAM, NIF#  NGP/THREE/A (3 ,5 1) , TOLD  (51)  ,TNEW(51) 

DO  2  1=1 , NGP 
2  TOLD (I) =TNEW (I) 

IF(MI.GT.O)  GO  TO  5 
11=1 

A  (2,  1)  =  1  .  +  ALPHA  (1)  *(1. +HTCKDZ(1) /COND  (1)1 

TOLD  ( 1)  =TOLD  ( 1)  ♦  ALPHA  ( 1)  *DZ  { 1 )  /COND  ( 1 ) * (HTC 1 *TF 1 +HS 1) 

GO  TO  7 
5  11  =  2 

TNEH  (1)  =32. 

TOLD  (II)  =TOLD  (I1)-A(1,I1)*32. 

7  IF(K.EQ.O)  GO  TO  12 
11=1 

DO  10  1=1, NLAM 
II=II*NDZ  (I) 

TOLD  (II)  =TOLD  (II)  ♦HRH  (1+  1)  *ALPHA  (I)  *DZ  (I)  /COND  (I) 

10  CONTINUE 

12  IF (M2.GT.0)  GO  TO  15 
I2=NGP 

A (2, NGP) =  1 .  ♦ ALPHA (NLAM) * ( 1 . *HTC2*DZ (NLAM) /COND ( NLAM) ) 

TOLD  (NGP)  =TOLD  (NGP)  ♦  ALPHA  (NLAM)  *DZ  (NLAM)  /COND  (NLAM)  *  (HTC:'*TF2+HS2) 
RETURN 
15  I2=NGP- 1 

TNEW  (NGP)  =  32. 

TOLD (12) =TOLD (12) -A(3,I2) *32. 

RETURN 

END 
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SUBROUTINE  GE06  (11,12,  A,  X,  t) 

C  SOLVES  A  SET  OP  THBEE-TEBM  LINEAR  EQUATIONS  USING  'GAUSSIAN 
C  ELIMINATION  AND  BACK  SUBSTITUTION 

C  REFERENCE:  APPLIED  NUMERICAL  METHODS  -  CARNAHAN  ET  ALIAE 
C  MODULE  NAME=GPS0« 

DIMENSION  A (3 ,5 1) , X (51)  ,  Y  (51)  ,  BETA  (51)  , GAMMA  (51) 

BETA  (II)  =A  (2,11) 

GAMMA  (I1)=Y  (II) /BETA  (II) 

111=11*1 
DO  2  1=111,12 

BETA  (I)  =A  (  2, 1)  -A  ( 1 , 1)  *  A  ( 3 , 1- 1)  /BETA  (I-  1) 

2  GAMMA  (I)  =(Y  (I) -A  (1,1)  *GAMMA  (1-1)  )  /BETA  (I) 

X (12) =GAMMA (12) 

1=12-11 
DO  5  K=1,L 
I=I2-K 

5  X(I)  =GAMMA  (I)  -A  (3,1)  *X  (I*1)/BETA  (I) 

RETURN 

END 


